%% Math 532 Project
% Andrew Price, Chris Millan, Crystal Bynum,
% Don Person, James Yerby, Spencer Carter

%% Housekeeping
clc;
clear;
pubOpts.evalCode = false;
publish('addEdge.m', pubOpts);
publish('alphaNormalize.m', pubOpts);
publish('CalculateGould.m', pubOpts);
publish('chooseBakeries.m', pubOpts);
publish('chooseBanks.m', pubOpts);
publish('chooseFireStations.m', pubOpts);
publish('chooseHospitals.m', pubOpts);
publish('choosePoliceStations.m', pubOpts);
publish('Coalesce.m', pubOpts);
publish('ColorEdges.m', pubOpts);
publish('ColorVertices.m', pubOpts);
publish('countSpanningTrees.m', pubOpts);
publish('covariance.m', pubOpts);
publish('Draw.m', pubOpts);
publish('DrawDG.m', pubOpts);
publish('GetDistances.m', pubOpts);
publish('GetEccentricities.m', pubOpts);
publish('GetMultipleLocations.m', pubOpts);
publish('GtoGstar.m', pubOpts);
publish('helpGraph.m', pubOpts);
publish('HighlightEdges.m', pubOpts);
publish('HighlightVertices.m', pubOpts);
publish('hurtGraph.m', pubOpts);
publish('isHamiltonian.m', pubOpts);
publish('isHamiltonianHelper.m', pubOpts);
publish('LineGraph.m', pubOpts);
publish('lineGraphDirector.m', pubOpts);
publish('makeHtmlTable.m', pubOpts);
publish('ParameterColor.m', pubOpts);
publish('randomGraphs.m', pubOpts);
publish('removeEdge.m', pubOpts);
publish('removeVertex.m', pubOpts);
publish('ShortestPath.m', pubOpts);
publish('toIncidence.m', pubOpts);
publish('toJIndexTable.m', pubOpts);
publish('wanderingDistanceMonteCarlo.m', pubOpts);


%% Load Matrices
load DigraphA.mat;
load K33Minor.mat;
%% Problem 1
%A
DrawDG(DigraphA);
title('Digraph A');
%B
%cmap = gray;
%cmap(:,1) = 1;
%makeHtmlTable(DigraphA,[],[],[],cmap);
spy(DigraphA);
title('Digraph A Adjacency Matrix');
%C
DigraphAIncidence = toIncidence(DigraphA);
figure;
spy(DigraphAIncidence);
%makeHtmlTable(DigraphAIncidence,[],[],[],cmap);
%% Problem 2
SCDG=zeros(60);
for i=1:60
    SCDG=DigraphA^i+SCDG;
end
SCDG;
%As can be seen from SCDG, which is the matrix DigraphA+DigraphA^2+
%DigraphA^3+...+A^60, some SCDG(i,j) = 0.
%Therefore, DG is not strongly connected.
%% Problem 3
%A
GraphA = cast((DigraphA | DigraphA'),'double');
Draw(GraphA)
%B
GraphA;
%C
%The incidence graph of DG is the same as the incidence graph as G.
%% Problem 4
SCG=zeros(60);
for i=1:60
    SCG=GraphA^i+SCG;
end
SCG;
%As can be seen from SCG, which is the matrix GraphA+GraphA^2+GraphA^3+...
%+GraphA^60, all SCDG(i,j) /= 0.
%Therefore, G is strongly connected. By strongly connected it is meant that
%you can reach any vertex of GraphA starting at any other vertex.
%% Problem 5
CoGraphA = Coalesce(GraphA);
%ids = ['41','31','18','39','30','20','42','40','11','6','21','26','27','36','15'];
bgG = biograph(tril(GraphA),[],'ShowArrows','off');
ids = [41 31 18 39 30 20 42 40 11 6 21 26 27 36 37];

bgG = (HighlightVertices(ids, bgG));

ids = {'Node 41' 'Node 31' 'Node 18' 'Node 39' 'Node 30' 'Node 20' ...
    'Node 42' 'Node 40' 'Node 11' 'Node 6' 'Node 21' 'Node 26' ...
    'Node 27' 'Node 36' 'Node 37' 'Node 39'};

edges = getedgesbynodeid(bgG, ids(1:15), ids(2:16));
% for i=1:max(size(ids))-1
%     bgG = ColorEdges(getedgesbynodeid(bgG, bgG.Nodes(i), bgG.Nodes(i+1)));
% end
set(edges, 'LineColor', [1 0 0]);
view(bgG);

view(biograph(tril(ones(15)) & K33Minor, ids(1:15),'ShowArrows','off'))
%G is not planar.
%The graph shown above is homeomorphic to K3,3 and is therefore not planar
%due to the Kuratowski Theorem stating that if a graph contains the
%subgraph K5 or K3,3 or a homeomorph of either it is not planar.

%% Computation of G*
%Pendant vertices- Vertices with degree 1.
%List- 1,12,13,19,24,34,35,38,43,47,49,52,55,56
GStar = GtoGstar(GraphA);
Draw(GStar);
%% Problem 6
sum(GStar);
%A Eulerian cycle exists IFF all vertices have an even degree and all
%vertices with a nonzero degree belong to a single contected component. G*
%has 14 vertices (4,9,10,17,18,25,28,29,30,31,37,46,48,50) with an odd 
%power. Thus add 7 edges connecting the odd powered vertices in pairs will 
%result in an Euler cycle making the graph Eulerian.
%% Problem 7
%Leave this segment commented to prevent long run times
%First we run the hamiltonian test on GStar.
%isHamiltonian(GStar)
%
%We note that GStar is not hamiltonian due to it not being biconnected
%After observing the graph of GStar, we can see that adding the following
%6 edges will create a hamiltonian cycle:
%22,23; 7,39; 11,16; 36,21; 45,42; 42,9
%G = GraphA;
%G = addEdge(22,23,G);
%G = addEdge(7,39,G);
%G = addEdge(11,16,G);
%G = addEdge(36,21,G);
%G = addEdge(45,42,G);
%G = addEdge(42,9,G);
%GStar = GtoGstar(G);
%
% Now we check to see if GStar is hamiltonian.  After a long exhaustive
% search, it is discovered to be hamiltonian
%isHamiltonian(GStar)

%% Problem 8
H2 = sort(sum(GraphA),'descend');
for i=1:60
    if H2(1,1) == 0
        %i-1
        break
    end
    for j=1:H2(1,1)
        H2(1,j+1) = H2(1,j+1)-1;
    end
    H2 = H2(2:60-i+1);
    H2 = sort(H2,'descend');
    H2;
end
%20 zeros
%% Problem 9
[LineG,LineGKey] = LineGraph(GraphA);
temp = LineGKey';
% Use the randomized L(G) generated by lineGraphDirector.m
% LineG = lineGraphDirector(LineG);
load directedLG.mat;
LineG = directedLineGraph;

for i=1:max(size(LineGKey))
    ids2{1,i} = strcat(num2str(temp(1,i)),',',num2str(temp(2,i)));     
end

%Number of edges in L(G)
sum(sum(tril(ones(max(size(LineG)))) & LineG));

%Draw Line Graph of G with correct labels
%DrawDG(LineG);
view(biograph(LineG,ids2))

%Incidence Matrix of L(G)
%C
LineGIncidence = toIncidence(LineG);

%% Problem 10
%For G, our range is 0.0602 >= p >= 0.0439
%For L(G), our range is 0.0443 >= p >= 0.0376
pg= randomGraphs(GraphA,0);
pl= randomGraphs(LineG,1);

%% Problem 11
distances = GetDistances(GraphA);
lgDistances = GetDistances(LineG);

%% Problem 12
[eccentricities, graphDiameter, graphRadius, ... 
    graphCenters, graphExtremes] = GetEccentricities(distances);

bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(graphCenters,  bgG, [1 0 0]);
bgG = ColorVertices(graphExtremes, bgG, [0 0 1]);

[lgEccentricities, lgDiameter, lgRadius, ... 
    lgCenters, lgExtremes] = GetEccentricities(lgDistances);

bgG = ColorEdges(lgCenters,  bgG, [1 0 0]);
bgG = ColorEdges(lgExtremes, bgG, [0 0 1]);
view(bgG)

%% Problem 13
[eccentricitiesBar, graphDiameterBar, graphRadiusBar, ... 
    graphCentersBar, graphExtremesBar] = GetEccentricities(distances, 1);

bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(graphCentersBar,  bgG, [1 0 0]);
bgG = ColorVertices(graphExtremesBar, bgG, [0 0 1]);

[lgEccentricitiesBar, lgDiameterBar, lgRadiusBar, ... 
    lgCentersBar, lgExtremesBar] = GetEccentricities(lgDistances, 1);

bgG = ColorEdges(lgCentersBar,  bgG, [1 0 0]);
bgG = ColorEdges(lgExtremesBar, bgG, [0 0 1]);
view(bgG)

%% Problem 14
%These functions will never actually be called due to time constraints
%since they have already been computed.

%wanderingDistance = wanderingDistanceMonteCarlo(GraphA,25000);
%lgWanderingDistance = wanderingDistanceMonteCarlo(LineG,25000);
load wanderingDistanceGMonteCarlo.mat;
load wanderingDistanceLGMonteCarlo.mat;
%Now for their eccentricities-
[wanderingEccentricities, wanderingGraphDiameter, wanderingGraphRadius, ... 
    wanderingGraphCenters, wanderingGraphExtremes] = GetEccentricities(wanderingDistance);
[lgWanderingEccentricities, lgWanderingDiameter, lgWanderingRadius, ... 
    lgWanderingCenters, lgWanderingExtremes] = GetEccentricities(lgWanderingDistance);

bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(wanderingGraphCenters,  bgG, [1 0 0]);
bgG = ColorVertices(wanderingGraphExtremes, bgG, [0 0 1]);

bgG = ColorEdges(lgWanderingCenters,  bgG, [1 0 0]);
bgG = ColorEdges(lgWanderingExtremes, bgG, [0 0 1]);
view(bgG)

[wanderingEccentricitiesBar, wanderingGraphDiameterBar, wanderingGraphRadiusBar, ... 
    wanderingGraphCentersBar, wanderingGraphExtremesBar] = GetEccentricities(wanderingDistance, 1);
[lgWanderingEccentricitiesBar, lgWanderingDiameterBar, lgWanderingRadiusBar, ... 
    lgWanderingCentersBar, lgWanderingExtremesBar] = GetEccentricities(lgWanderingDistance, 1);


%% Problem 15
%Here we will normalize each of the specified indeces (e(x),Id(x),we(x),
%and Iwd(x)) with the a(x) -> a'(x) as specified in the problem
%calculations for G
eccentricitiesPrime = alphaNormalize(eccentricities);
eccentricitiesBarPrime = alphaNormalize(eccentricitiesBar);
wanderingEccentricitiesPrime = alphaNormalize(wanderingEccentricities);
wanderingEccentricitiesBarPrime = alphaNormalize(wanderingEccentricitiesBar);
%calculations for L(G)
lgEccentricitiesPrime = alphaNormalize(lgEccentricities);
lgEccentricitiesBarPrime = alphaNormalize(lgEccentricitiesBar);
lgWanderingEccentricitiesPrime = alphaNormalize(lgWanderingEccentricities);
lgWanderingEccentricitiesBarPrime = alphaNormalize(lgWanderingEccentricitiesBar);

%% Problem 16
%Compute the variances of the different eccentricity variables
%Calculations for G
varEccentricitiesPrime = var(eccentricitiesPrime);
varEccentricitiesBarPrime = var(eccentricitiesBarPrime);
varWanderingEccentricitiesPrime = var(wanderingEccentricitiesPrime);
varWanderingEccentricitiesBarPrime = var(wanderingEccentricitiesBarPrime);
%Calculations for L(G)
varlgEccentricitiesPrime = var(lgEccentricitiesPrime);
varlgEccentricitiesBarPrime = var(lgEccentricitiesBarPrime);
varlgWanderingEccentricitiesPrime = var(lgWanderingEccentricitiesPrime);
varlgWanderingEccentricitiesBarPrime = var(lgWanderingEccentricitiesBarPrime);
%Compute pairwise covariances for each combination
%Calculations for G
covGeId = covariance(eccentricitiesPrime,eccentricitiesBarPrime);
covGewe = covariance(eccentricitiesPrime,wanderingEccentricitiesPrime);
covGeIwd = covariance(eccentricitiesPrime,wanderingEccentricitiesBarPrime);
covGIdwe = covariance(eccentricitiesBarPrime,wanderingEccentricitiesPrime);
covGIdIwd = covariance(eccentricitiesBarPrime,wanderingEccentricitiesBarPrime);
covGweIwd = covariance(wanderingEccentricitiesPrime,wanderingEccentricitiesBarPrime);

CovarTable = [1 covGeId covGewe covGeIwd;
              covGeId 1 covGIdwe covGIdIwd;
              covGewe covGIdwe 1 covGweIwd;
              covGeIwd covGIdIwd covGweIwd 1];
makeHtmlTable(CovarTable);

%Calculations for L(G)
covLGeId = covariance(lgEccentricitiesPrime,lgEccentricitiesBarPrime);
covLGewe = covariance(lgEccentricitiesPrime,lgWanderingEccentricitiesPrime);
covLGeIwd = covariance(lgEccentricitiesPrime,lgWanderingEccentricitiesBarPrime);
covLGIdwe = covariance(lgEccentricitiesBarPrime,lgWanderingEccentricitiesPrime);
covLGIdIwd = covariance(lgEccentricitiesBarPrime,lgWanderingEccentricitiesBarPrime);
covLGweIwd = covariance(lgWanderingEccentricitiesPrime,lgWanderingEccentricitiesBarPrime);

%% Problem 17
%The values for a,b,c,d need to be changed in toJIndexTable.m
%G
jGIndex = toJIndexTable(eccentricitiesPrime,eccentricitiesBarPrime,wanderingEccentricitiesPrime,wanderingEccentricitiesBarPrime);
jGIndexPrime = alphaNormalize(jGIndex);

bgG = biograph(tril(GraphA),[],'ShowArrows','off');
ParameterColor(eccentricities, bgG);
%ParameterColor(eccentricitiesBar, bgG);
ParameterColor(wanderingEccentricities, bgG);
%ParameterColor(wanderingEccentricitiesBar, bgG);
ParameterColor(jGIndex, bgG);
%L(G)
jLGIndex = toJIndexTable(lgEccentricitiesPrime,lgEccentricitiesBarPrime,lgWanderingEccentricitiesPrime,lgWanderingEccentricitiesBarPrime);
jLGIndexPrime = alphaNormalize(jLGIndex);

%% Problem 18
%G
jGDiameter = max(jGIndexPrime);
jGRadius = min(jGIndexPrime);
jGCenters = (jGIndexPrime == jGRadius);
jGExtremes = (jGIndexPrime == jGDiameter);

bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(jGCenters,  bgG, [1 0 0]);
bgG = ColorVertices(jGExtremes, bgG, [0 0 1]);

%L(G)
jLGDiameter = max(jLGIndexPrime);
jLGRadius = min(jLGIndexPrime);
jLGCenters = (jLGIndexPrime == jLGRadius);
jLGExtremes = (jLGIndexPrime == jLGDiameter);

bgG = ColorEdges(jLGCenters,  bgG, [1 0 0]);
bgG = ColorEdges(jLGExtremes, bgG, [0 0 1]);
view(bgG)

%% Problem 19
%[V,D]=eig(GraphA)
%v1=V.^2;
GouldIndices = CalculateGould(GraphA);
GouldIndicesPrime = alphaNormalize(abs(GouldIndices));
GouldjGCompare(1,:) = jGIndexPrime;
GouldjGCompare(2,:) = GouldIndicesPrime;
makeHtmlTable(GouldjGCompare);

%% Problem 22
disp('Bakery Locations');
bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bakeries = chooseBakeries(distances, wanderingDistance, 3);
bgG = HighlightVertices(bakeries(1), bgG);

%% Problem 23
bgG = ColorVertices(bakeries(2), bgG, [0 1 0]);
bgG = ColorVertices(bakeries(3), bgG, [0 0 1]);
view(bgG);

%% Problem 24

%disp('Bank Locations');
banks = chooseBanks(distances, wanderingDistance, 3, 2);
bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(banks(1), bgG);
bgG = ColorVertices(banks(2), bgG, [0 1 0]);
bgG = ColorVertices(banks(3), bgG, [0 0 1]);
view(bgG);

%disp('Fire Station Locations');
fireStations = chooseFireStations(distances, 3);
bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(fireStations(1), bgG);
bgG = ColorVertices(fireStations(2), bgG, [0 1 0]);
bgG = ColorVertices(fireStations(3), bgG, [0 0 1]);
view(bgG);

%disp('Police Station Locations');
policeStations = choosePoliceStations(distances, wanderingEccentricities, 5);
bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(policeStations(1), bgG);
bgG = ColorVertices(policeStations(2), bgG, [0 1 0]);
bgG = ColorVertices(policeStations(3), bgG, [0 0 1]);
view(bgG);

%disp('Hospital Locations');
bgG = biograph(tril(GraphA),[],'ShowArrows','off');
bgG = ColorVertices(GetMultipleLocations(distances, 3), bgG, [0 0 1]);
bgG = ColorVertices(GetMultipleLocations(distances, 2), bgG, [0 1 0]);
bgG = ColorVertices(GetMultipleLocations(distances, 1), bgG);
view(bgG);

%% Problem 25
%publish('DrawLineGraphs.m');
DrawLineGraphs;

%% Problem 26
%By using Kirchhoff's matrix tree theorem, we find that L(G) has
%6.5879e+036 spanning trees, which is significantly greater than
%Avogadro's number (6.0221415 e+023)
numSpanningTrees = countSpanningTrees(LineG);
